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Abstract. Theories of monochromatic high-frequency electromagnetic fields 
have been designed by Felsen, Kravtsov, Ludwig and others with a view to por- 
traying features that are ignored by geometrical optics. These theories have 
recourse to eikonals that encode information on both phase and amplitude — 
in other words, are complex- valued. The following mathematical principle is 
ultimately behind the scenes: any geometric optical eikonal, which conven- 
tional rays engender in some light region, can be consistently continued in the 
shadow region beyond the relevant caustic, provided an alternative eikonal, 
endowed with a non-zero imaginary part, comes on stage. 

In the present paper we explore such a principle in dimension 2. We in- 
vestigate a partial differential system that governs the real and the imaginary 
parts of complex- valued two-dimensional eikonals, and an initial value problem 
germane to it. In physical terms, the problem in hand amounts to detecting 
waves that rise beside, but on the dark side of, a given caustic. In mathe- 
matical terms, such a problem shows two main peculiarities: on the one hand, 
degeneracy near the initial curve; on the other hand, ill-posedness in the sense 
of Hadamard. We benefit from using a number of technical devices: hodograph 
transforms, artificial viscosity, and a suitable discretization. Approximate dif- 
ferentiation and a parody of the quasi-reversibility method are also involved. 
We offer an algorithm that restrains instability and produces effective approx- 
imate solutions. 



1. Introduction 

1.1. Geometrical optics fits well a variety of issues, but especially survives as 
an asymptotic theory of monochromatic high-frequency electromagnetic fields - 

[BB] . [Du] . [FKN] . [US] . [33, Hi], |KE] > [EH] and [K32], (KPT] . [EH], [HQ], [E3 

are selected apropos references. Generalizations have been worked out by Felsen, 
Kravtsov, Ludwig and their followers — see e.g. [CF1] and [CF2] , [EFj . [Fel| and 
[Fi2], [HF], [KrT] . [Kr2] . [KTH] and [KFA] . [LBL] . [LuT] and [Lu2], or consult [BM] . 
[CLOT] . [KQ2j . One is enough for successfully modeling basic optical processes, 
such as the propagation of light and the development of caustics. The others 
embrace geometrical optics and are additionally apt to account for certain optical 
phenomena — for instance, the rise of evanescent waves past a caustic — that are 
beyond the reach of geometrical optics. A leitmotif of these is allowing a keynote 
parameter to adjust itself to a standard equation, and simultaneously take complex 
values. 
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(1.1) 



The following partial differential equation 

dw 
dy 



( dw 
\~dx~ 



n 2 (x,y) 



underlies the mentioned theories in case the spacial dimension is 2. Here x and y 
denote rectangular coordinates in the Euclidean plane; n is a real-valued, strictly 
positive function of x and y; w is allowed to take both real and complex values. Func- 
tion n represents the refractive index of an appropriate (isotropic, non-conducting) 
two-dimensional medium — its reciprocal stands for velocity of propagation. Func- 
tion w is named eikonal according to usage, and relates to the asymptotic behavior 
of an electromagnetic field as the wave number grows large — the real part of 
w accounts for oscillations, the imaginary part of w accounts for damping. Geo- 
metrical optics deals exclusively with real- valued eikonals, complex- valued eikonals 
distinguish a more advanced context. Throughout the present paper we assume the 
refractive index is conveniently smooth, and consider sufficiently smooth eikonals. 
The following partial differential system 



(1.2) 



= n 2 (x,y) 



^x^x 



UyVy = 



governs those complex-valued solutions to (jl.ip that obey 

u = Rew, v = Imro. 

Observe the architecture of (II. 2|) : gradients are involved through their orthogonal 
invariants — lengths and inner product — only. Observe also the following prop- 
erties, which result from a standard test and easy algebraic manipulations. First, 
system (|1.2p is elliptic-parabolic or degenerate elliptic. Second, a solution array [u v] 
to (|1.2p is elliptic if and only if its latter entry v is free from critical points. 
The Backlund transformation, which relates u and v thus 



(1.3) 



and implies both 



and 



r.2 



-1 

1 



Vu, 



|Vu| : 



sgn / = sgn (u x v y - u y v x ), 



\Vu\ > n 



(1.4) 




Vu } = 0, 



is another, decoupled form of JT7 

System (|1.2[) discloses two scenarios — the former is tantamount to conventional 
geometrical optics, the latter opens up new vistas. Either the following equations 



2,2 2 

u x + u y = n 



and 



= Vy=0 



hold, or the following inequalities 

|Vu| > n and |Vd > 0, 
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and the following equations prevail 
Vu= 1 J 



(1.5) 



1 

-1 



Vv, 4 = 1 




(1.7) (|Vu| 4 — n 2 u 2 )u xx + 2n 2 u x u y u X y + (|Vu| 4 — n 2 u 2 )u yy — n\Vu\ 2 (Vn, Vu), 

(1.8) (|Vv| 4 + n 2 ^)^ - 2n 2 v x v y v xy + (| Vw| 4 + n 2 v 2 x )v yy + n\ Vw| 2 (Vn, Vu) = 0. 

Equations (| 1 . 5|) invert the Backlund transformation mentioned above, and imply 
(jl.6p . (jl.7[) and (|1.8p are quasi-linear partial differential equations of the second 
order in non-divergence form. The former has a mixed elliptic-hyperbolic character: 
a solution u is elliptic or hyperbolic depending on whether the length of V« exceeds, 
or is smaller than n. The latter is elliptic-parabolic or degenerate elliptic: a solution 
v such that Vw is free from zeros is elliptic, degeneracy occurs at the critical points 
of v. 

Any real- valued sufficiently smooth solution to (| 1 .4|) satisfies (jl.7l) . A real- valued 
function is an elliptic solution to (11. 7|) if and only if it coincides with the former 
entry of an elliptic solution to (|1.2|) . A real- valued function u, smooth and without 
critical points, is a hyperbolic solution to (|1.7p if and only if two real-valued smooth 
functions ip and ip exist such that 

and 

u = -(ip + tp). 

Any real- valued, sufficiently smooth solution to (| 1 . 6[) satisfies (|1.8p . Any eZ%>- 
t?c solution to (II. 8|) satisfies (|1.6|) . A solution to (|1.8|) need not satisfy ljl.6|) : for 
instance, perfectly smooth solutions to (|1.8|) exist whose gradient vanishes exclu- 
sively in a set of measure 0, and which make the left-hand side of (|1.6p a non-zero 
distribution. 

Let J be endowed with an appropriate domain and obey 

J(v) = //i(^) 

for any u from that domain — here j is the arc length along a parabola, videlicet 

j(p) = |Vi + p 2 + 2 lo §(p + Vi + p 2 ) 

for any real p. The following properties hold, (i) J is convex, coercive and sub- 
differentiable, but not Frechet-differentiable. (ii) Any critical point of J, i.e. any 
function v such that 

dJ(v) 3 0, 

satisfies (|1.6p in any open set O such that 

O is essentially contained in {(x, y) S domain of v : Vt>(a;, y) 7^ 0}. 
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Consequently, a critical point of J solves a free-boundary problem for equation (| 1 .6|) 
— the relevant free boundary is 

(domain of v) PI d{(x, y) G domain of v : Vv(x, y) =/= 0} 

and plays the role of a caustic, (ii) Any function v such that 

J(v) = minimum 



satisfies (jl.8p in an appropriate viscosity sense. In other words, a minimizcr of J 
solves in a generalized sense a boundary value problem for equation (JTT8J) . 

An early treatment of (|1 .2|) traces back to [ER]. Further apropos information 
is offered in [MTT] , [MT2] . [MT3] , and |MT4j . where solutions in closed form, 
qualitative features, exterior boundary value problems, related free boundaries, 
variational and viscosity methods are discussed. 

1.2. Geometrical optics ultimately amounts to manipulating: (i) the Riemannian 
metric known as travel time, videlicet 



n(x, y)y/dx 2 + dy 2 ; 

(ii) the members of appropriate one-parameter families of travel time geodesies — 
called rays; (iii) the envelopes of rays — called caustics. 

Geometric optical eikonals are entirely controlled by rays. They shine in light 
regions (those spanned by relevant rays), burn out beside caustics (where the ray 
system breaks down), and shut down in shadow regions (that rays avoid). As a 
consequence, geometrical optics is unable to account for any optical process that 
takes place beyond a caustic, on the dark side of it. Essentials of two-dimensional 
geometrical optics, which are instrumental here, are outlined in Appendix A. 

Complex- valued eikonals are more flexible. As the cited work of Felsen, Kravtsov 
and Ludwig may prompt, complex- valued eikonals look apt to consistently continue 
geometric optical cognates into shadow regions. Such a continuation is the main 
theme of the present paper. 

Let us pave the way by heuristically considering the case where refractive index n 
is 1. A classical recipe informs how general geometric optical eikonals can be cooked 
up: start from a complete integral, derive a one-parameter family of solutions, take 
the relevant envelope, and shake well. Let / be an arbitrary, but sufficiently smooth, 
real function. The following pair 

(1.9) w = x cost + y sin t + fit), = — xsint + y cost + fit) 

causes w and t to enjoy the following properties: 



dw \ 2 / dw N 2 



dx ) \ dy ' 

the pertinent eikonal equation governing w; 

dt . dt n 
cost — — h suit — =0, 
ox ay 

a Burgers-type equation governing t; 

W X t X + Wyty = 0, 

showing that the gradients of w and t are orthogonal; 

w x — cost, w v — sint, 
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a Bdcklund transformation further relating w and t. Both the rays of w and the 
level lines of t are the straight dines where 

— xsint + ycost + f'(t) = 

and t equals a constant. Such straight-lines span the light region and envelope the 
caustic. We have 

— /' = the support function of the caustic, 

and 

x = /"(t) cost + /'(£) sin t, y = f"(t) sin t - f'{t) cost, w = f(t) + f"(t) 

along the caustic. Therefore the second-order derivatives of w and the gradient of t 
simultaneously blow up there — in particular, the caustic of w is also the shockdine 
oft. 

We claim that both w and t can be continued beyond the caustic, in a subset of 
the shadow region, if complex values are allowed. Suppose / is analytic — so as it 
can be continued by a holomorphic function of a complex variable. Let i = \J — 1, 
the unit imaginary number. Let u, v, A, /i be real; put the following equations 

w = u + iv, t = A + i/_t, 

and equations Q1.9P together, but force Q1.9P to produce real x and y. The following 
formulas 

sin A , , . cos A , , , cos A , , . sin A , , . 

x = — — Re/' (t) + — — Im/' (t , y = — Re/' (t + — — Im/' t), 

cosn /x sinn /i cosh /x sinn /z 

u = Re/'(t) + (a; cos A + ysin A) cosh /x, v = Im/'(t) + (—a; sin A + ycos A) sinh/x, 

ensue, then an inspection testifies that the claimed continuation ensues too. 

Incidentally, we have also shown that solutions to the non- viscous Burgers equa- 
tion can be continued past the shock-line by suitable complex-valued solutions to 
the same equation. Let us mention that complex-valued solutions to viscous and 
non- viscous Burgers equation are dealt with in |K WYZ| . |PSj and [KOj . A more 
exhaustive analysis is carried out in Appendix B. 



1.3. A certain initial value problem for system (II. 2[) — the one described in items 
(i) and (ii) below, and in figure 1 — summarizes our purposes. 

(i) An initial curve IC is given. The following alternative applies: either IC is 
specified exactly — no error infects the definition of IC; or else IC is a phantom — 
some coarse and polluted sampling of IC is gotten. In the former case assume IC 
is smooth enough. In the latter case recover IC, i.e. fed the available data into an 
appropriate denoising process, and then elect the consequent output as an operative 
substitute of IC — ad hoc tools are provided in Appendix C. 

Represent IC (either the authentic one, or else its surrogate) by the following 
equations 

(1.10) x = a(t), y = P(t), 

and adjust parameter t so as 

t = a travel time 



without any loss of generality. 
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Figure 1. A geometric optical eikonal and its continuation past a caustic. 



Assume travel time is an extra metric in action and the relevant geodesic curva- 
ture of IC is free from zeroes. In other words, postulate that (| 1 . 1 0|) and either of 
the following equations 

k (velocity) 3 = Geodesic curvature, 
k (velocity) 2 = Euclidean curvature — (unit normal, V logn(a;, y)) 

result in 

k vanishes nowhere. 

(ii) A pair [uv] is sought that obeys system f| 1 . 2[) and fulfills the following conditions. 
First, 

(1.11) u(x,y)=t, v(x,y)=0, 

if x, y and t are subject to (|1.10|) . Second, u and v are defined in the side of IC that 

(sgn k) x (unit normal to IC) 

points to. 

Arguments from Appendix A allow us to comment as follows. IC and the men- 
tioned side of it can be viewed as a caustic and a shadow region, respectively. Any 
geometric optical eikonal, which makes IC a caustic, lives in the illuminated side 
of IC; the complex-valued eikonal, whose real part is u and whose imaginary part 
is v, lives in the opposite, dark side of IC; both equal a travel time along IC. An 
extension of the geometric optical eikonal in hand ensues. Such an extension does 
obey the eikonal equation, lives in both the light region and a subset of the shadow 
region, and takes complex values where shadow prevails. In physical terms, prob- 
lem (i) and (ii) accepts a caustic in input, then models evanescent waves that rise 
in the dark side of it. 
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In the present paper we focus our attention on solutions [u v] to the problem (i) 
and (ii) that meet the following extra requirements: 

(iii) they are elliptic; 

(iv) their Jacobian determinant obeys 

Sgll (U x Vy ~ UyV X ) = SgRK. 

Condition (iii) ensures that the latter entry v is not constant; as will emerge from 
subsequent developments, initial conditions plus conditions (iii) and (iv) ensure 
that the same entry is nonnegative. 

The solutions to problem (i)-(iv) develop singularities near IC, as any geometric 
optical eikonal does in the vicinity of the relevant caustic. We will show that they 
obey the following expansions 

(1.12) u(x, y) = s + o(r), v(x, y) = ^|«( a )|V3| r |3/a + o(r 3/ 2) 

as (x, y) belongs to the appropriate side of, and is close enough to IC. Here r and 
s stand for the curvilinear coordinates described in Appendix A — informally, r 
is a signed distance from IC, s is a lifting of a travel time inherent to IC. Note 
the physical meaning of (jl . 12[) : the damping effects, which are encoded in the 
imaginary part of a complex- valued eikonal, are tuned by the geodesic curvature of 
the relevant caustic. 

2. Framework 

2.1. A convenient coordinate system must be called for. We choose to recast (|1.2[) 
by reversing the roles of dependent and independent variables — i.e. we think of 
u and v as curvilinear coordinates, and think of x and y as functions of u and v. 
In other words, we subject (|1.2p to the change of variables that is called hodograph 
transformation at times — see [Evj and Zw , for instance. 

Let [u v] be any smooth elliptic solution to (|1.2p , and observe the following. 

(i) The level lines of u and those of v are free from singular points, and cross at a 
right angle. Moreover, the Jacobian determinant 

(2.1) u x v y — u v v x vanishes nowhere. 
Let 

[uv] i — * [x(u, v) y(u, v)] be a local inverse of [xy] i— > [u(x, y) v(x, y)], 
then observe the following equation 

{u x v y - u y v x ){x u y v - x v y u ) = 1 

and the propositions (ii)-(iv) below. 

(ii) The following partial differential system holds 

(2.2) 1/E-1/G=l, F = 0. 
Here 

E = n 2 (x,y)(xl + y 2 u ), F = n 2 (x, y)(x u x v + y u y v ), G = n 2 (x, y){x 2 v + y 2 L ,) 
— in other words, 

n 2 (x,y)((dx) 2 + {dy) 2 ) = E{du) 2 + 2Fdudv + G{dv) 2 . 
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(iii) The following systems and equations hold 

d r x 

du y 



(2.3) 
(2.4) 



/ 



d_ 

dv 



x 
v 






1 " 


d 


X 


-1 





dv 


. y . 


" 


-1 " 


d 


X 


1 





du 





(2.5) 



sgn/ = sgn(x u y v - x v y u ), 



f 2 



1 + n 2 ( Xl y){x 2 v + y 2 v ), f 2 = l- n 2 (x, y)(x 2 u + y 2 u ). 



(iv) The following equations hold 

1 - n 2 (x,y)xl 



(2.6) 



d 2 



r 



1 + n 2 (x, y)x 2 



d 2 
dv 2 



x 

y 



= (W 2 ) 



1 -n 2 (x,y)yl 
1 + n 2 (x,y)y 2 ' 

x u ~~ Vu ~~2< x uyii 

2a; u y u x u — y u 



V\ogn(x,y) 



(2.8) 



in the event that n is identically 1, these equations read 
1 — x 



? \ 2 ?\2 



J %w ~ 0; Uuu ~f~ 



i + vl 



Vvv = 0. 



l + x 2 

Proof of (i). System (11.21) tells us that 

|V«| > 0, 

and that the gradients of u and v are orthogonal. An assumption gives 

|Vu| > 0. 

The following equation 

(u x v y - u y v x ) 2 = |Vu| 2 |Vv| 2 - (Vu, Vv) 2 
concludes the proof. □ 



Proof of (ii). Since 



d(u, v) 



d(u,v) d(x,y) 
d(x,y) d(u,v) 



d(x,y) 
d(x,y) 



d(u, v) 



we have 



1 



|V«| 



<9(it, v) 


T 




_d(x,y)_ 






d(x,y) 






d(u, v) 




X 


(Vu, 


Wv} 2 



1 
1 

|Vu| 2 (Vu.Vu) 
(Vu,Vv) \Vv\ 2 



x u ~l~ y u x u x v -\- yuyv 

^ _1_ «, », -r-2 I „,2 



1 



|v,| 2 (Vw ' Vw> 



^ 2 + vl |Vw| 2 ' a;2+y2 |y u |2 

XuXv +y u y v = -(Vu,Vv)(x u y v - x v y u ) 2 . 
The conclusion ensues. □ 

Proof of (iii) and (iv). The latter equation from (|2.2[) yields 

Xu • yv — yu • ( x v)i 
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hence the following systems result 

Xy, p 1 

Vu ~ J -io 



Xy 



Xy 

Vv 



X u 

Vu 



Factor / is easily identified. Both the above systems give 

sgn/ = sgn(x u yv - x v y u ); 
the former and system (|2.2[) imply 

-p = l + n 2 (x,y)(xl+yl); 

the latter and system (|2.2j) imply 

f 2 = l-n 2 (x,y)(x 2 u +y 2 u ). 
Another arrangement reads 










/ " 




Vu 




Vu 




-f~ 






Xy 




. -v/ 







Vv 


1 


Vv 




. i// o 




Xy 



— two Backlund transformations, inverse of one another. The former and system 
([23 imply 

1 - n 2 (x,y)yl 



r = 

the latter and system (|2.2|) imply 



1 



1 



n 2 {x,y)yl' 

n 2 (x,y)xl 



1 



n 2 {x,y)x% ' 

The integrability conditions, which pertain to the Backlund transformations in 
hand, read 

—] \ l /f 1 [ Xu Vu ] = 

du dv / 

I x v y v j 

and result in equation (|2.7|l after algebraic manipulations. 

Equations (|2.5p and (|2.6[) are consistent with one another and with the early 
equations (|1.3p and (| 1 . 5[> . as Proposition (ii) and its proof show. 

The proof is complete. □ 

2.2. In view of the discussion above, problem (i)-(iv) stated in Subsection 1.3 
amounts to looking for solutions [x y] to the following partial differential system 

d_ 

(2.9) dv 
sgn/ = sgn k, f 2 = l-n 2 (x,y)(x 2 u +yl), 

that are defined either in the half-strip 

a < u < b, 0<t><oo 
or in an appropriate bounded subset if it, and satisfy the following initial conditions 

(2.10) x{u,0)=a(u), y(u, 0) = 0(u) for a < u < b. 

We will concentrate on such a problem. A behavior of solutions [x y] to (|2.9p and 
(|2.10| as v is close to is fixed up in the next section. An algorithm for computing 
the same solutions is offered in Section 4. Section 5 is devoted to an example. 



X 


i 


' -1 " 


d 


X 


. y . 


= 7 


1 


du 
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3. Behavior near the caustic 



The state of affairs causes any solution of (|2.9j) and (|2.10p to suffer from singu- 
larities near the initial line — indeed, 

x 2 v {u,v) +yl(u,v) -> oo 

as a < u < b and v approaches 0. The proposition below offers more details on the 
subject, as well as a proof of expansions (|f .f 21) . 

Proposition 3.1. Let x and y obey system (|2.9[) and initial conditions (|2 . 10|) . As- 
sume x(u, v) and y(u, v) depend smoothly upon u for every nonnegative, sufficiently 
small v. Then the following asymptotic expansion holds 



x(u, v) 




a(u) 


_ y(u,v) 




/?(«) 



(3u) 2 / 3 sgn n(u) 



-f3'(u) 
a!{u) 



+ o(v 2 / 3 ) 



2\ K (u)\ 1 / 3 ^a'{u) 2 + f3'(u) 2 
as a < u < b and v is positive and approaches 0. 

Proof. A hypothesis made on a and (3 in Subsection f .3, equations (|2.9p and initial 
conditions (|2.10|) tell us that 

/(«,«) ^0 

as a < u < b and v is positive and approaches 0. Therefore, 
„i/3 



lim ■ 



Inn sgn flu. r) <jsgn/(u,u) / -^-f 3 (u,v) 



vlO f{u,v) vlO 

by L'Hospital's rule. Equations (|2.9j) give successively 

7r n ( x >y) = — Tt — r 

ov f[u,v) 
d 



1/3 



V log n(x,y), 



dv 



K + vi) 



f{u,v) 



and 



dv 



f(u,v)=3(xl + y 2 u )V 2 n 2 (x,y) 



{xl + 2/2)3/2 
Vlogn(x,y),(x 2 u +yl)- 1/2 



We infer 



lim (3v) 1/3 



Vv 



sgn k(u) 



because of the very definition of k. 
The conclusion follows. 

4. Discrete setting 



-f3'(u) 
a'(u) 



a 



4.L Rendering problem ()2.9|) and (12. f 0j> into an effective discrete form entails 
coping with singularities of solutions, overflows, and ill-posedness in the sense of 
Hadamard. 

Singularities result from features of both the system and the initial conditions in 
hand, as already remarked in the previous section. Overflows take place whenever 
the constraint 

n 2 (x,y){x 2 u + y 2 u ) < I 
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chances to be violated. Note that the system 



d 


X 




1 " 


d 


X 


dv 


. y _ 




-1 


du 


y 



- Cauchy-Ricmann, a possible linearized version of (|2.9 
solutions 

x 

y 

which are highly instable, i.e. 



possesses obvious 



= exp(— vt + tv) 



cos(iu) 
sin (tit) 



sup ■ 



-i 2 



+ 



-oo < u < oo 







-oo < u < oo 



} 



as t ] oo and j = 1, 2, . . . , and 

inf {x 2 (u, v) + y 2 (u, v) 

as t j oo and v is positive. 

Ill-posedness is a typical drawback of initial value problems for partial differential 
equations and systems of elliptic type. It was observed by Hadamard [Halj - Ha2j. 
and deeply investigated by John [Jol] - |Jo3] . Lavrentiev jr. |Lvlj - |Lv3j . [LVj . Miller 
[Mll]-[Mi2], Payne [FaI]-[Fa4| and [PSTJ IESg, Pucci [Pcl]-[Pc4], and Tikhonov 
[Til] - [Ti5] . Classical surveys on the subject have been authored by Lavrentiev jr., 
[Lv4] . [LRS] , Payne [Pa5] . and Tikhonov [H]. Related information is in [ErT]-|Er2]. 
[BLW] . [BG] . [E], [TVT] . [Mr], [Ha], [E] and [ ^. A sample of more recent contri- 
butions includes [M], [Al], [SKI, [578], pCP], [BE] , [BR] . [BV] . [E], [Ci], 
El], [CD], [CE], [C3] , [Cd] , [Pi] . [EUl-lEll], [EwT]-[Ew2]. [HH] . [HR], [Hg, [KV] . 
[Kn] . [LV] . El, [MaT]-[Ma2], [Mo] . [Ntl]-[Nt2], [Pi6]-[Pa8], [EE], [RE], [RHH] . [RS] . 
[Rm], [St], [TGSYj . [TLY] . [VFC] . [Waj . 

As experience suggests, one might attempt to contend with ill-posedness via a 
priori bounds on solutions and similar devices. We opt not to touch on this issue 
in the present paper, and focus our attention on constructive aspects instead. 

We rely upon: (i) asymptotic expansions, describing how relevant solutions be- 
have near the initial line; (ii) a technique of approximate differentiation, especially 
designed for working in presence of errors; (iii) an appropriate injection of arti- 
ficial viscosity, which softens a coefficient and protects against overflows; (iv) an 
imitation of the quasi-reversibility method. 

4.2. Besides data exposed to view, our method takes six parameters in input: 
M, N, A, [i, v, £. The first two are large integers that specify the number of samples 
in hand — e.g. M — N = 100. The remaining four parameters set the tone of 
smoothing processes: A and fi relate to approximate differentiation; v stands for 
viscosity; £ relates to quasi-reversibility. 

4.3. The following equations and inequality 

ustep = (b — a)/{M — 1), vstep > 0, 
uj = a + (j — 1) ustep (j = 1, . . . , M), Vk = (k — 1) vstep (k — 1, . . . , N), 

will be in force throughout this section. We choose a mesh to consist of the following 
points 



(uj,v k ) (j = 1, 



,M:k = 1, 



,N), 
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and store sample values at mesh points in the following matrices 

[ X (J> fe)]i=l,...,M;fc=l,...,iV, [V(j> k )]j=i,...,M-k=l,...,N- 

The columns of these matrices, i.e. 

x(;l),x(;2),...,x(;N), y(; 1), y(; 2), . . . , y(; N) 

are recursively generated in the following way. 
First step. 

1) = at(uj), y(j, 1) = P(uj) (i = 1, . . • , M), 

according to initial conditions (|2.10p . 
Second step. 

sgnn( Uj ) ~P'(uj) 



5(7,2) =x(j,l) + 



(3v 2 ) 



2/3 



y{j, 2) - y(j, 1) + 



sgn n{v,j) a'juj) 

2|«(u i )|V3 ^a'( % -) 2 +/3'(%0 ; 



(3« 2 ) 2/3 (j = l,...,M), 



according to expansion (|3.ip . 

Further steps. For fc = 3, ... ,7V do actions (i)-(iii) below, 
(i) Differentiate. 

X = X (;k-l), Y = y(;k-1), 

A = DX, B = DY. 



Here 



A = a dimensionless positive parameter, 

— — ~ 0.0415 + 0.5416 x M ~ 1/2 - 0.6426 x M" 1 + 1.3706 x M~ 3/2 , 
b — a 



K{u) 



1 



cos{ut) 



dt 



for every u, and 



D = - 
t 1 



K 



/' 



Aid 



i,i=l,...,M 



A' 



/' 



i,j=l,...,M, 



— a matrix that mimics differentiation with respect to variable u, and is analyzed 
in Appendix C. 
(ii) Enter viscosity. 

/(j) = P( V ,n(X(j),Y(j)) y/A{j)* + sgn/s^) (j = 1,...,M), 



£7 



Here 



/(I) o 
/(2) 





5, 1/ 



/(i) o 

/(2) 




•• f(M) 
< i/ = artificial viscosity < tt/2, 








A. 
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and P is given by 



sinv 



cos^ v 



if < p < (1 +sin 2 z/T 1 / 2 , 
if p> (l + sin 2 ^)- 1 ^ 



— observe that 

< p P(y, p) is strictly positive and continuously differentiable, 
P(v, p) approaches \f \ — p 2 uniformly as < p < 1 and v approaches 0. 
(iii) Enter quasi-reversibility. 

x(-,k) = ip(v k ), y(-,k) = ip(v k ). 
Here ip and ip are the vector-valued mappings that are generated thus: 
£ = a dimensionless positive parameter; 



R = (ustep)' 



2 -5 4 -1 ••• 

1-2 1 ••• 

••• 1-2 1 

••• -1 4 -5 2 



a caricature of a second-order derivative; 

<p{vk-i) = X, <p'(v k -i) = U, 



3t\\R^v k )\\ 2 + (v S tep) J y"(v)\\ 2 dv 



minimum, 



VK-i) = r, ^'K-i) = v, 

ft 

3£ ||i?V(vfe)|| 2 + (wsiep) y ||V>"0)|| 2 cfo = minimum. 

Vk-l 

As is easy to see, 

tp""(v) = if w fe _i < w < v k , 
iff\v k ) = 0, ^'"K) = H{vstep)- 2 {R T R)p{v k ); 

4)""(v) = if < v < v k , 

i>"(v k ) = 0, tp'"(v k ) = 3Z(vstep)- 2 (R T R)i>(v k ). 

Therefore, 

x(-,k) = (Id + (,vstep(R T R)y 1 (X + vstepU), 

y(-,k) = (Id + £ vstep {R T R))^ 1 (Y + vstepV). 

in other words x(-,k) and y(-,k) are mollified versions of X + vstepU and 
y + vstepV, respectively. 
Last step. End. 
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4.4. As a matter of fact, the above process simulates the following partial differ- 
ential system 



d 


X 


i 


' 


-i " 


d ' 


dv 


. y _ 


~7 


i 





du 



<9u 4 



/ = P (u, n(x, y)yjxl + y£j sgnn(u), 



where the modified, and extra, terms protect against overflows and instability. 
The methods based on either artificial viscosity or quasi-reversibility share basic 
features: they all suggest perturbing the underlying partial differential equation or 
system in a way or another, in order to palliate obstructions. The quasi-reversibility 
method was introduced in |LL] , and improved in Mi5j . |GZ| ; other references are 
[AT] . [AtY] . [BiT]-[B32], [DR], [EwT] . po] . [rTu] . [HZ] . [KS] . [Li], [Pil], |5EI| - 
[Sh2] . [TT] . 



5. Example 

For simplicity, suppose refractive index n is 1. Consider the curve — known 
as Tschirnhausen's cubic or trisectrix of Catalan |Lwj . |Se] — whose parametric 
equations read 

x = ^(1 — 3t 2 ), y=^(3 — £ 2 ), (t — parameter), 



and imply 



arc length = — (t 2 + 3), t = 2sinh ( -arcsinh(arc length) 



oo 4 ctooo * 



* o 

o 



3oo 08«ooooO°' 



-3 -2 
x-axis 



Figure 2. An arc of Tschirnhausen's cubic: a gross sampling 
(stars) and a denoised version (cicles). 




x-axis 



Figure 3. Level lines where either u = constant or v = constant. 



A geometric optical eikonal w making Tschirnhausen's cubic a caustic is repre- 
sented by the following equations 

1, 9N 2st t, 9N s(l-t 2 ) 
x=-(l-3t 2 ) y= -(3-t 2 ) + - 

w = s + —(3 + t 2 ) (s, t — parameters) 

in the light region. As arguments from Appendix B show, the same eikonal can be 
continued in the shadow region via the following equations 

1 - 2(,s 2 + i 2 ) + s 4 - 2s 2 t 2 - 3i 4 * [3 - 2(s 2 - t 2 ) - (s 2 + t 2 ) 2 } 



2(1 + s 2 +t 2 ) ' V ~ 2(1 + s 2 + t 2 ) 

2(is + t) l-(is + t) 2 

l + (is + t) 2X+ l + (is + t) 2 



W = is + 1 — - — — — — % + - | / . — ' y (s,t = real parameters). 



The method from Sections 1 to 4 goes in the following way. 
(i) Consider the arc of the Tschirnhausen's cubic where 

-1.5 < t < 2.2, 

and let 

(j l 3D 

be a gross sampling of such an arc — in other words, 

i i = -1.5 + 3.7^ (j = l,...,31), 
Uj = i(l - it 2 ) + (5% random noise), 0j — ^-(3 — t 2 ) + (5% random noise). 
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Figure 4. Plot of imaginary part v versus x and y. 



(ii) Plug gross data into the following denoising process 

A = 0.005, n = 0.1260, 



31 

E 

3=1 
31 

E 

3=1 



.7-1 



/3 



30 



+ A / [//(a"") 2 +^ 1 a 2 } dt = minimum, 



+ A / [//(/? ) +fi~ L p z ]dt = minimum, 



and let the path that is represented by the following equations 

x = a(t), y = /3(t), < i < 1, 

surrogate the original Tschirnhausen's cubic, 
(iii) Adjust parametric equations as follows: 

dt 



du 



a'(t) 2 + p'(t) 2 }- 1/2 , *(0)=0, 



x = a(t(u)), y = (3(t(u)), < u < Length. 

(iv) Select requisite parameters thus 

M = 101, N = 91, vstep = 0.005, v = 0.5, £_ = 0.9; 

and then set the algorithm from Section 4 to work. 

Results are shown in figures 2, 3 and 4, and comfortably agree with those drawn 
from closed formulas. 
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Appendix A. 

Basic mathematical lineaments of two-dimensional geometrical optics are out- 
lined in the next paragraphs. Selected references on geometrical optics, and on 
some of its generalizations and applications are |BB| . |Duj . |GS| . [Jnj . [Kej . |KL| . 

[Em-IKI^, IKKl, IE3H, P^], IHEI, 

A.l. Terminology. Let n be a refractive index — i.e. a tractable function of two 
real variables x and y, which takes positive values only and is bounded away from 
zero locally. Any real-valued, suitably smooth solution w to is a geometric 

optical eikonal (GOE). The domain of w, plus parts of the relevant boundary, is 
a light region; the complement of it is a shadow region. The trajectories of Vw, 
namely the orbits of the following dynamical system 

dx dy 
w x (x,y) w y (x,y) 

are called lines of steepest descent — irrespective of whether they are genuine lines 
or not. A line of steepest descent is a ray if w is twice continuously differentiable 
in some neighborhood of it; any line of steepest descent, which is not a ray, is a 
caustic. (Rays are smooth curves, which have one degree of freedom and travel all 
over areas without intersecting one another. Caustics are exceptions in a sense: 
loosely speaking, they can be thought of as envelopes of rays.) The Riemannian 
arc length, whose element takes the following form 

n(x, y) \J dx 2 + dy 2 , 

is known as travel time — travel time is an alias of the customary arc length in the 
case where n = 1. 

A. 2. GOEs and travel time, (i) The restriction of any GOE w to either an 
appertaining ray or caustic automatically coincides with a properly rescaled travel 
time t. 

(ii) Let nodal line be an alias of locus of zeros. The value of any GOE w at any 
point (x, y) equals either the travel time between (x, y) and a nodal line of w or the 
negative of such a travel time — provided (x, y) is not a long way off. 
Proof of (i). By definition, both rays and caustics of w obey 

dx : w x (x, y) = dy : w y {x, y); 

the following equation 

t = a rescaled travel time 

is an alias of 

dt = ±n(x, y)^J(dx) 2 + (dy) 2 . 

Consequently, 

(dw/dt) 2 =n- 2 (x,y)(w 2 x + w 2 y ) 

along any ray or caustic in question. Property (i) follows. □ 

Property (ii) is a consequence of (i) and Fermat's principle below. 
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A. 3. Fermat's principle. The travel time geodesies, i.e. those curves which ren- 
der 



J n{x,y) yjdx 2 + dy 2 



either a minimum or stationary, are characterized by the following second-order 
ordinary differential equation 

curvature = (unit normal, V log n(x, y)) 

— they have geodesic curvature 0, and are perfect straight lines in the case where 
n = 1. The rays of any GOE are geodesies with respect to travel time. The travel 
time geodesies that are trajectories of some proper vector field — i.e. have one 
degree of freedom and are free from mutual intersections — are the rays of some 
GOE. 

The foregoing statements rest upon first principles of calculus of variations, dif- 
ferential geometry and ordinary differential equations. Recall that the following 
formulas apply to any smooth parametric curve: 



unit tangent = (velocity) 
-id 



velocity = \J{dx/dt) 2 + (dy/dt) 
i d 



dt 



unit normal 



II -f 
1 



(unit tangent) 



(velocity) —(unit tangent) = (curvature) x (unit normal), 



(velocity) 1 — (unit normal) 

curvature = (velocity) -3 

curvature = (velocity) -3 /). -3 



-(curvature) x (unit tangent), 

dx/dt d 2 x/dt 2 
dy/dt d 2 y/dt 



x(t - h) y(t - h) 1 
x(t) y(t) 1 

x(t + h) y(t + h) 1 



+ 0(h 2 ), 



curvature = (velocity) 3 h 



x(t) y(t) 1 

x(t + h) y(t + h) 1 
x(t + 2h) y(t + 2h) 1 

x(t + h) y(t + h) 1 
x(t + 2h) y{t + 2h) 1 
x(t + 3h) y(t + 3h) 1 



0(h 2 ). 



Recall that the following formula 

curvature of the lines of steepest descent = 
Jacobian determinant of |Vw| -1 & w 

applies whenever w is smooth and has no critical point. Recall also that 

n(x,y) x (Geodesic curvature) = 
Euclidean curvature — (unit normal, V log n(x, y)) 
if travel time is an alternative Ricmannian metric in force. 
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A. 4. Initial value problems and geometry of their solutions. The condition 
of taking given values along some given path is qualified initial according to usage. 
Seeking a GOE, which obeys some initial condition, is an initial value problem. 
Such a problem has either two different solutions or no solution at all, depending on 
whether the eikonal equation and the initial condition match or conflict. Generally 
speaking, a solution w can be detected in the former case by successively detecting 
the objects listed below, based upon the arguments provided. 

• The values of Vw along the initial curve. 

Since the eikonal equation specifies the length of Vw and the initial condition 
identifies a tangential component of Vw, the normal component of Vw along the 
initial curve comes out in two different modes. 

• The rays of w. 

An ODE reasoning demonstrates that the travel time geodesies, which live near the 
initial curve and leave it with the same direction as Vw, are the trajectories of a 
smooth vector field. By Format's principle, these geodesies are the requested rays 
indeed. 

• The values of w itself on each ray. 
Property (i) from Subsection 2 fits the situation well. 

A. 5. Standard initial value problems. The present item concerns existence, 
regularity and the number of GOEs that satisfy orthodox initial conditions. 

Assume henceforth all ingredients are smooth and let IC stand for initial curve 
in shorthand. Let 

(A.l) x = a(t), y = P(t), a<t<b 

be a parametric representation of IC such that 
(A. 2) t = a travel time. 

Let the initial condition imply 

w{x,y) = 7(i) 
as x, y and t satisfy (jA.lj) . and assume 



(A.3) 



< 1 



for a < t < b. Then exactly two GOEs satisfy the initial condition displayed. 
Moreover, these eikonals are smooth in a full neighborhood of IC — the relevant 
light regions surround it completely. 

The case where refractive index n is constant involves explicit formulas, as well as 
gives evidence to interactions among the eikonal equation, Burgers-type equations 
and Backlund transformations. Assume 

n = 1, 

and let (IA. 1[) to (|A.3|) be in force. Define ip, ip,ui,p, q by 

cost/? = a', sin <p = f3'; 



cos4> — j', sintfj = ±yT — (7O 2 ; 
to = tp + tp: 
p = cosw, q = smui. 
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Since the appurtenant Jacobian matrix equals 



q(s) p(s) 
-p(s) q(s) 



sin^(s) — (w — 7(s)) duj(s)/ds 




the following pair 



X 




a(s) 


. y _ 




I3(s) 



+ (to - 7(a)) 



p{s) 



makes s and w implicit functions of x and y in a neighborhood of IC. The properties 
listed below ensue. Function s obeys the following equations 

-q(s) {x - a{s)) +p(s) (y - p(s)) = 0, 



p(s) 



5s 
dx 



ds 
dy 



0. 



dx ) 



( 2i 

dy 



smip(s) — (w — 7(s)) 



duj(s) 



ds 



(The first assures that the level lines of s are straight, the second is a PDE of 
Burgers type.) Function w obeys the following eikonal equation 



dw 
dx 



dw 
dy 



1. 



Functions s and w are related by the following equations 



^xx ^xy 


dto(s) 


-Wy 


V^xy ^yy 


ds 


w x 



(The former is a Bdcklund transformation, which pairs solutions to the Burgers and 
the eikonal equations mentioned above. It assures that Vs and Vro are orthogonal, 
and the level lines of s are both lines of steepest descent and isoclines of w.) The 
Euclidean metric obeys 

dx 2 + dy 2 = |Vsp 2 ds 2 + dw 2 . 

There holds 

s(x,y) = t, w(x,y) = j(t) 
as x,y and t satisfy (|A. 1[) . If uj is free from critical points, the line where 

smtp(s) 



X 




a(s) 


. y _ 




. /3( fi ) . 



diu(s)/ds 



p(s) 



is a caustic. (Such a line is an envelope of the level lines of s. Function s develops 
shocks along it, the restriction of w to it equals arc length, and the second-order 
derivatives of w blow up there.) 
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A. 6. Borderline initial value problems and caustics. The present item is a 
recipe for producing caustics, which involves coupling the eikonal equation with 
borderline initial conditions. In the case where refractive index n is identically 1, 
any smooth convex curve can be viewed as a caustic provided a GOE is detected, 
whose restriction to the curve in hand equals a relevant arc length. 

Let (|A.1|) specify IC and let (|A.2|) hold. Let the initial condition imply 

w(x,y) = t 

as x, y and t satisfy (|A.1|) . Assume that the appropriate geodesic curvature of IC is 
free from zeros — in other words, let (|A.1|) and the following equation 

k (velocity) 2 = Euclidean curvature — (unit normal, V logn(a;, y)) 



result in 



K vanishes nowhere. 



Then exactly two GOEs satisfy the present initial condition. Both these eikonals 
fail to exist on both sides of, and be smooth near IC. They turn IC into a caustic, 
and make the side of it, which 

(sgn k) x (unit normal) 

points to, a shadow region. Either eikonal in hand obeys 



2 v / 2 3 i 
w(x,y) = s± — — |r| = |k(s)| 5 + 0(r 2 



Aw(x,y) = ± -^=\ r \ 



1*001* + 0(1) 



at every point (x, y) that belongs to the light region and is close enough to IC — in 
particular, the two-sheeted surface, made up of the two eikonals in hand, exhibits 
an edge of regression above IC. 

Here r and s are the curvilinear coordinates that the following pair 



X 




a(s) 


. y _ 







+ r(a'(s) 2 + (3'(sf 



a'(s) 



relates to rectilinear coordinates x and y. Coordinate r is a signed distance from 
IC, coordinate s makes (a(s(x, y)), f3(s(x, y))) the orthogonal projection of (x,y) 
on IC. The former is constant on the parallel lines to IC and obeys the following 
eikonal equation 

2 / a \ 2 

= 1; 



dr 
dx 



the latter is constant on the normal straight-lines to IC and obeys the following 
Burgers-type equation 

ds/dx ds/dy 
a'(s) P'{s) 

Both are subject to the following Backlund transformation 



= 0. 



Vr = [a'{sf +(5'{sf 



-P'(s) 
a'(s) 



and exhibit singularities along the evolute of IC. 
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Appendix B. 

Here we sketch a method of continuing a two-dimensional geometric optical 
eikonal past a caustic. Though rigorous, such a method is slightly reminiscent 
of the so-called theory of complex rays — cf. [CLOT] or [KFA] . for instance. It 
applies in the case where the refractive index equals 1, and can be used for si- 
multaneously continuing solutions to non-viscous Burgers equation beyond shock 
lines. 



B.l. Our method involves analytic continuation from the real- number axis into 
the complex plane — an ill-posed process in the sense of Hadamard. Let h be a 
real or complex-valued function of a real variable, or even a list of samples. An 
analytic continuation of h is a holomorphic function of a complex variable, whose 
domain surrounds the real axis and whose restriction to the real axis fits h well — in 
other words, a solution H of the following initial value problem for Cauchy-Ricmann 
equation 

— , H(.,0)ch. 

If h is an analytic function, and is not polluted by noise, an analytic continuation 
H of h results from obvious formulas. For example, 

k=0 



oo 

H{x, y) = l- J exp[i£(ai + iy)] h{£) d£, 



where hat denotes Fourier transformation. 

If h collects gross data, an effective analytic continuation H of h can be obtained 
by analytically continuing an appropriate, smoothed and denoised version of h. 
Consider for instance the case where 

— oo < a < b < oo, 
N = an integer larger than 1, 

a mesh size is given thus 

Ax = (b — a)/(N - 1), 

nodes are given by 

Xj = a+(j-l)Ax {j = 1,...,N); 

and h stands for 

hi, ■ ■ ■ , hN, 

a sequence of noisy samples. 

An ad hoc analytic continuation H solves the following least square problem 

N b+Ax/2 L 

Ax ^^[H(xj ,0) — hj] 2 + -^-jr J dx J \H(x, y)\ 2 dy = minimum 

■ J = 1 a-Ax/2 -L 



2:S 



in a convenient class of holomorphic functions — e.g. the class of trigonometric 
polynomials of a suitable degree. Here 

b+Ax/2 L 



A J dx J \H(x,y)\ 2 dy 



a-Ax/2 —L 

plays the role of a penalty; A and L are regulating parameters — A is related to 
noise, L is related to a priori information. 

An explicit expression of H can be derived via discrete Fourier transforms. Sup- 
pose for simplicity that N is odd, say 

iV = 2n + l (n = 1,2,...); 

let 

T = N Ax, 

and let DFT be the discrete Fourier transform that obeys the following equations 

N 

DFT k {h) = hj exp (-2mk^ (k = -n, . . .,n), 
1 n 

h i=N £ DFT »( h ) ex P ( 27 ™^) (i = 1, • • • , N), 

k——n 

n N 



k——n j—1 



cf. IBHI. for instance. If 



then 



sinh(47rfc£/T) 
= ' fe = 47r fcx/r ' C - k =Ck ( fc = !' • ■ ■ ' n )> 

H(x, y)= + XC k )- l DFT k {h) exp Umk^^\ . 



k- 

This is a T-periodic trigonometric polynomial of degree n that enjoys the following 
properties 

N / , x 2 N 



and 



El^o)-M 2 <( XTI7 ^J EN 2 > 

b+Ai/2 L 

± J dx J \H(x,y)fdy<^, 

a-Ax/2 —L 

N / _i_ * \ 



— here Dn denotes the Dirichlet, or periodic sine function obeying 

sin(A^/2) 
Dn{x) = ATsin(x/2) 

if x/(2ir) is not an integer. 
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More information on analytic continuation can be found in |Abj . |A12j . |BVlj . 

[EE], H [DSU-M. & EES], Em], [S3, Ei, 

[Evgj . [LAI - [Li], [LRS], [Mi3]-[Ml], [MY], El, E3, [Sy], PS], US], [Ve], [S3, 

EH. 



B.2. Consider a plane curve C that either is inherently smooth or results from a 
suitable smoothing process of raw data. Assume C is analytic and its curvature 
vanishes nowhere. For simplicity, assume C is the graph of the following equation 

V = f(x), 

and / is convex. 

Alternative parametric representations of C, which are instrumental throughout, 
include 

(B.l) x=t, y = f(t), 

where parameter t coincides with the abscissa; and 
(B.2) x = g'(t), y = tg'(t)-g(t), 

where parameter t is the slope of the tangent straight-line. Here g denotes the 
Legendre conjugate of / — recall from e.g. [Roj that / and g are related thus 

t = f(x), x = g'(t), tx = f(x)+g(t), 1 = f"(x) g"{t). 

Curve C changes into a caustic under the following modus operandi. Let 

(B.3) x = a(t), y = (3(t) 

be any parametric representation of C, where a and (3 are analytic. Let 7 and k 
stand for arc length and curvature, respectively — in other words, 



7 (i) = / Va'(t) 2 + 0'(t) 2 dt, 
a"/3' 



a' 13" 
[(a') 2 - 



The following pair 
(B.4) 



X 




a(s) 


. y . 




f3(s) 



-(/3')2]3/2- 

w — 7(s) 



a'(s) 



y/a'{s) 2 +f3'{s) 2 

makes s and w curvilinear coordinates. As is easy to see, the level lines of s are 
tangent straight-lines to C, the level lines of w are involutes of C — orthogonal to 
one another. 

Equations (|B.4p imply 

y < fix), 

and 

s(x,y)=t, w(x,y)=j(t), 

ifx,y and t obey (|B.3|) — s and w live below C and satisfy precise conditions along 
C. 

We compute 



d{x,y) 
d(s, w) 



-f3'{s) a'(s) 
a'(s) p'(s) 



k(s)[w — 7(s)] 

" [(a') 2 + (f3')T 1/2 
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to draw the following set: 
B.5) 

B.6) 

B.7) 

B.8) 
B.9) 
B.10) 



dw\ 2 ^ ( dw\ 2 
dx] \dy J 

ds dw ds dw 



1, 



0, 



dx dx dy dy 
Vw = [a'(s) 2 +0'(s) 2 ]- 1 / 2 
\Ws\- 2 = [w- 1 ( S )} 2 K(s) 2 [a'(s) 2 +p'(s) 2 }, 



a'(s) 
(3'(s) 



W X y 



1 



7(s) 



-w v 



-Wy W X \ 



Equation (|B.5|) is a conservation law; it reads 



ds 



0, 



9s 

the standard Burgers equation, if (|B.2p is in force. Equation (|B.6[) is the equation 
of geometrical optics in hand. Equation (|B.7|) shows that the gradients of s and w 
are orthogonal. Equation (|B.8|) shows that both C and the tangent straight-lines to 
C are lines of steepest descent of w; it also shows that the straight-lines in question 
are isoclines of w. Equation (|B.8j) can be viewed as a Backlund transformation, 
which converts any solution to (|B.5[) into a solution to (|B.6[) . It reads 

1 

/'(*) 



v w = [i + /'( s ) 2 ]- 1 / 2 



s = 9 — 

w x 



or simply 



-1/2 



1 



depending on whether (|B.1|) or (|B.2|) is in effect. Equations (|B.9p and (|B.10[) show 
that both the gradient of s and the second-order derivatives of w blow up near C. 

We infer that s is governed by a Burgers-type equation, and develops shocks 
along C. The following objects — w, C, the tangent straight-lines to C, and the 
region below C — are a geometric optical eikonal, the relevant caustic, the rays, 
and the light region, respectively. 

We now claim: (i) s and w can be continuously extended into the region where 

(B.ll) y > f(x), 

the dark side of C, if suitable imaginary parts are provided; (ii) the relevant exten- 
sions obey equations (|B.5[) to (|B.10[) . 

The points above C are reached by no tangent straight-line to C, of course. We 
insist in drawing tangent straight-lines from these points, but allow complex slopes. 
In other words, we recast (|B.4[) this way 

-f3'(s)[x-a(s)]+a'(s)[y-l3(s)]=0, 
[w - 7(*)]K(s) 2 + P'is) 2 ] 1 ' 2 = a'(s)[x - a(s)] + f3'(s)[y - /3(a)], 
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and force such equations to hold in the situation where 

Re(s) = A, Im(s) = //, Im(a;) = Im(y) = 0. 
The following formulas result 



_ Im[a>{s)(a{s)[3'(s)-a'(s)l3(s))} _Im\J3'{s)(j3(s)a'(s)-(3'(s)a(s))] 
x _ ~~ > y — 



Jm[a'(s)P'(s)] 



Jm{j3'(s)a'(s)} 



(B.12) 

w = 7(3) — 



s = X + i/A, 

lmf3(s) 



a'(s) 



Im[a'(s)/3'(s)] 



■13' (a) 



Ima(s) 



lm[(3'(s)a'(s)} 



where a, j3 and 7 stand for the analytic continuations of the original objects. 
Formulas (|B.12|) answer the claim. Among other things, they give 



f + o(m 4 ), 



X 






dA(X) 


V(A) " 


H 2 + B(X) 


-/3'(A) " 


. y _ 




. pw . 


dX 


J3'(X)_ 


_ a'(A) _ 



A 



1 log| K |-ilogV(a') 2 + (/5') 2 , 5 



(/?') 



3 01 1 2 
as /i approaches zero, and 

, d(x,y) \a'{s)Imf3(s)- f3'(s)Ima(s)\ 2 \a / (s)/3"(s)-a"(s)f3 / (s)\ 2 

det dix^) = ' 

consequently 
and 



Ima'(s)/3'(s) 



1/ - fix) = - /"(a(A)) [Im a(A + ^)] 2 + 0( M 4 ) 
det = /u [a'(A)/3"(A) - a"{X){3'{X) + 0{^ 2 )} 



9(A, M ) 

as /i approaches zero. Thus (|B.12p imply (|B.11[) . as well as 

d(X,n) 

if is different from, and sufficiently close to 0. 

B.3. Here is an example. A catenary is the graph of either the following equation 

y = cosh x 

or the following equations 

(B.13) x = log(t+ \Jl + t 2 ), y=^/l + t 2 , 

where parameter t coincides with both an arc length and the slope of the tangent 
straight-line. 

Consider solutions to the following equations 

dw\ 2 /<9uA 2 ^ ds ds 
dx J \dy J ' dx dy 
which obey the following conditions 

w(x,y) = s(x,y) = t 
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as x and y obey (IB.13j) . Function w and the catenary in question are an eikonal 
and the relevant caustic, respectively; s obeys Burgers equation, takes a constant 
value on each tangent straight-line to the catenary and develops shocks along the 
catenary. 

The light region is the set where 

— oo < x < oo, y < coshx; 

the shadow region is the set above the catenary, where 

— oo < x < oo, cosh a: < y. 

The following pair 



s x — y — s log(s + yl + s 2 ) — \fl + . 



w 



+ sy - log(s + vl+ s 2 ) 



VTTT 2 

makes s and w implicit functions of x and y in the light region. Eikonal w and 
its partner s can be continued in a subset of the shadow region via the following 
equations 

x = A + tanh A (fi cot fi — 1), y = — - [ sinh 2 A h /i sin /i + cos fi 

cosh A \ sin fj, 

s = sinh(A + ifx), 

• i ^ M i f ■ ^ 

w = smn A — 1 — - (sin /i — fi cos /i) 

sin cosh A 

— here A and fi are parameters such that 

-oo <A<oo, 0</i< ir/2. 
The foregoing equations imply that 

~1 



y — cosh x — (cosh A) ji 



2+0(M 2 ) 



as fi approaches 0. Furthermore, 

d(x, y) [fi 2 + tanh 2 A (1 — fi cot fi) 2 ) (sinh 2 A + cos 2 fi) 



dct 



d(X, fi) cosh A sin fi 

dw 1 dw , . . 

— — = tanh(A + z/i) 



9a; cosh(A + ifi) ' 9y 
— in particular, a singularity occurs at the point whose coordinates are 

x = 0, y = 7r/2. 
Figures 5 and 6 show plots of the imaginary parts of w and s. 

Appendix C. 

C.l. Differentiating a real- valued function of one real variable is among the most 
elementary processes of mathematical and numerical analysis, but is also a signifi- 
cant prototype of those problems that are nowadays qualified ill-posed in the sense 
of Hadamard. Methods of approximating derivatives of smooth functions under 
non-exact data have been widely experimented over the years. Here we take the 
opportunity of sketching one more of such methods. We consider the case where 
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Figure 6. Burgers equation: the imaginary part of s beyond a 
shock-line. 



data consist of discrete and noisy samples, nodes are equally spaced, and informa- 
tion is available on both the relevant noise and the underlying smoothness. Our 
method is inspired by ideas that the theory of statistical learning has recently re- 
vived — see e.g. [US], [EPP] . [SZl]-[SZ4] 
several of its ascendants — see e.g. |ACR] . 

[cjw], ECU, UM, EDo], pp. [Egg], [EU, m 

[LW] , [MM] . [MtiTj - fMul] . [Mg] , [mmz] . [qT 



Vp4l 



[S3, M, [Wn]- 



and of course mimics 
MT]- [XB2] . |5H] . [Ba] , [BrT] . 
3E], [KM], [KW], [Kq|, [LP] - 
[RS]-[Ra2], [RS], [Eg, [Sk], [SL], 
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Items in input include: 

(i) the end points of a bounded interval — a and 6; 

(ii) the number of both nodes and samples — an integer N, larger than 2; 

(iii) nodes from a to b — specifically, 

x k =a+(k-l)^± (k = l,...,N); 

(iv) noisy samples — a sequence 

of real numbers whatever. 

Goals include recovering some function / and the derivative /' of / based upon 
the following information only: 

(v) / is smooth; 

(vi) f{xk) is close to gk, for k = 1, 2, . . . , N. 

Our recipe segments into the following three steps. First, let A and p, be positive 
parameters, and solve the following variational problem 

N °9 
(C.l) y~" y [u(x k ) - gk] 2 + A / [pL 7 {u"") 2 + p^u 2 ] dx = minimum, 

under the following condition 

(C.2) u belongs to Sobolev space W 4 ' 2 (—oo, oo). 

Second, adjust A and fi properly. Third, take u, u' as approximations of /, /'. 

C.2. Effective formulas read as follows. 

Rudiments of the calculus of variations demonstrate that Problem (jC.lj) & ((C.2)) 
possesses a unique solution; moreover that such a solution — named u throughout 
— obeys 

( d 8 u \ N 

(C.3) xl^-^+^uj +J2l u ( x k)~ gk]5(x-x k ) = 

^ X ' k=l 

for — oo < x < oo. The following features are decisive: equation (IC3[) is affine; 

" dx^ +fi 

is a positive operator in L 2 (~ oo, oo), whose inverse mollifies; 

N 

^2[u(x k ) -gk]S(- -x k ) 
fc=i 

is a spifce iram or a shah-function. 

Condition (|C.2[) and equation (|C.3[) give 

iv 

(C.4) -Au(i) = ^[^(x fc ) - 5fc ]A- 1 ' ~ ' 1 

k=l 
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for — oo < x < oo. Equation (IC.4j) gives 
(C.5) A 



u{xi) 




u(xi) 




9i 




+ A 




= A 




u(x N ) 




_ u(x N ) 




. 9n _ 




Figure 7. Plots of K,K' and K". 

Here A denotes an appropriate fundamental solution to d 8 /dx s + 1, namely the 
solution to 



that decays at infinity — A is given by 



(C.6) 



7r A (a:) 



cos(x£) 



for — oo < a; < oo, is even and positive definite. Secondly, 



.4 



K 



K 



A(0) 



x 2 - x\ 



x N — X\ 



A 



K 



x\ - x 2 
A 4 

A(0) 

£jy - ^2 
A* 



A' 



A' 



A 4 

^2 - %N 

A* 

A(0) 



a symmetric, positive definite Toeplitz matrix. 
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Let v = 7r/8. Manipulating formula (|C.6j) gives 
1 la; 2 1 x 4 1 



jr(a?) 



1 \x\< 



isinz^ 



1 cos v 24 



! sin 720 2 5040 



+ 0(x 8 ) 



as x approaches zero — K behaves near zero like a spline of order seven. Calculus 
of residues, or convenient formulas from [PBM, Section 2.5.10], give 



4K(x) = e-W^cosi 



|x| sin j/ - i/) + e ^ smu sm(\x\ cosz/ + v) 



as — oo < x < oo. Therefore, 



d n K 



dx 



-(x) 



(— sgnx) r 



| e -Wc««' C0S (| a .| 8i]ai/ _ + 



e - |x|sin!/ sin(|x| cosi. + (n + l)v - wr/2)} 
as x ^ and n = 1, 2, 3, • • • . Figure 7 shows plots of K, K\ K" . 




10 12 14 16 18 20 



Figure 8. Reciprocal condition estimator of A, plotted versus the 
ratio (b-a)/((N -l)fi). 



Analysis shows the following. The spectrum of A lies in the open interval 
]0, iV/ (8 sin i/) [. All eigenvalues of A are close to l/(8sinz/), if fi(N — l)/(6 — a) 
is small; otherwise, the largest eigenvalue of A is close to AT/ (8 sin i/) and most 
remaining eigenvalues of A are close to 0. In particular, A is invertible anyway; A is 
either well- conditioned or ill-conditioned depending on whether fi(N — a) is 

small or large. Figure 8 shows plots of a reciprocal condition estimator of A versus 
(b-a)/((N -!)»). 

Equation (|C.4[) implies that u belongs to the linear span of 
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— translations and dilations of K. Such items own either a spike-shaped or a well- 
rounded profile depending on whether /i is small or large, inasmuch as 



SZo \ (d/dx)Item\ 2 dx _ 5 



— fi tan v, 



\(d 2 /dx 2 )Item\ 2 dx 
|Item| 2 dx 



— u tan v. 
7 



The same items are definitely linearly independent, although appropriate formulas 
and analysis show that their Gram matrix is well-conditioned only if fj,(N—l) / (b—a) 
is small enough. 

Equation (|C.5j) determines u(xi), . . . ,u(xn) in terms of data. Note that they 
solve 



N 



y^}u(x k ) - gk} 2 + A [u(xi) ■ ■ ■ u(x N )} A 1 



k=i 



u{xi) 
u(x N ) 



minimum 



— a standard finite- dimensional least-square problem, where A and the inverse of 
A imitate a regulating parameter a la Tikhonov and a penalty, respectively. 

Now we are in a position to draw conclusions. Let Id be the N x N unit matrix, 
and 

R = (Aid + Ay 1 
— a resolvent. Let B the vector-valued function such that 



B(x) = R 



K 



K 



x — X\ 



x — xn 



for —00 < x < 00 — an alternative basis in the linear span mentioned above. Define 
C and D thus 

C = AR 



D = il- 



ly' 



, ( x 2 - Xi 



The following equations hold. 



for —00 < x < 00, 



K 



, xi-x 2 



, I X N - Xi \ , ( X N - X 2 



= [9i ■■■ 9n] B(x) 



u(xi) 




9i 




= C 




_ u(x N ) 




. 9N _ 



K 



fJL 



K 



, x 2 -x N 



/' 



R. 
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u'(xi) 




9i 




= D 




_ u'{x N ) 




. 9n _ 



0.6 



0.5 - 



0.4 



0.3 




asse T 

Figure 9. The alternative basis B : plots of B\{x), . . . , B N (x) 
versus x. 




Figure 10. Matrix C playing the role of a regularizing filter. 

As figures 9, 10 and 11 show, B and C mimic a typical Green's function from 
two different perspectives, D mimics a derivative of a Green's function. Observe 
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Figure 11. Matrix D, simulating differentiation. 



incidentally that 



J [pL 7 (B"") T B"" + l i- 1 B T B]dx=tx [A (Aid + A)~ 2 ] , C = Id - A (Aid + A)' 1 , 

— oo 

and that B and C solve the following variational problem 



A 



:3C 

j [ti 7 (B"") T B"" + f i- 1 B T B] dx + tr [(C T - Id)(C - Id)] 



minimum, 



subject to the following conditions 

Be \W**(-oo,oo)] N , C = 



Bx( Xl ) B 2 { Xl ) 
Si (12) B 2 (x 2 ) 

Bi(xn) B 2 (xn) 



B n (xi) 
B N {x 2 ) 

B n (xn) 



C.3. Here we offer directions for adjusting A and /1 properly. 

Parameter A — dimensionless — discriminates whether solution u to Problem 
(|C.1[) & (|C.2p either fits data well (but is simultaneously sensible of noise), or else 
is little affected by noise (but departs somewhat from data). In loose terms, the 
following statements hold. First, u virtually interpolates g\,.. - ,gjv if A is close to 
zero; however, u is liable to own an irregular profile at the same time. Second, u 
quenches smoothly if A grows larger and larger. Indeed, 
00 

J [J(B"") T B"" + M - X S T B] dx = tr A- 1 + 0(A), 



C = Id + 0(A) 
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as A approaches zero; 



\d n B{x)/dx n \\ < X- 1 



N 



8sin[(n + l)v] 

as — oo < x < oo and n = 0, 1, . . . , 6. 

Parameter [i — making fi/(b — a) dimensionless — determines how much u is 
close to, or departs from a spike train. Indeed, 

5(x — x\) 



B(x) = /i • 



1 + A sin v 

as — oo < x < oo and fi approaches zero; 

1 



S(x — xn) 



B(x) = 



N + 8A sin v 



1 



as — oo < x < oo and /i grows larger and larger. 




-0.2 0.2 

logarithm of rho 



Figure 12. Plots of S.R.Q. versus log p. Here p = (6 - a)/((N - l)/i). 



In the case where suitable information is available a priori, a theorem from 
Subsection IC.4I below suggests which values of A and /i work properly. Other- 
wise, parameter A may be determined based upon the discrepancy principle, the 
cross-validation, the L-curve criterion, or other customary devices — see e.g. |Hn[ 
Chapter 7] for details. 

Parameter fi is expediently identified by the following recipe. Let 



R.Q. of v = /i 1 



^[(d/dxYvjx^dx 
.Coo Hx)] 2 dx 
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a dimensionless Rayleigh quotient; define a relevant minimum thus 



S.R.Q. = min < R.Q. of v : ^ v e span of K 



K 



then determine /i so that 

(C.7) S.R.Q. = minimum. 

Equation (|C.7[) causes solution u of Problem (|C.1[) & (|C2[) to retain its most 
favorable Rayleigh quotient. Equation (|C.7[) can be approached via equation (|C.8|) 
below and tools from linear algebra, although care must be taken of the ill-condition 
of the involved matrices. 

Let 



G = /i 1 x Gram matrix of K 
H = fi^ 1 x Gram matrix of 



■ - Xi 
A* 

• - xi 



,K 



x N 



\ l J V A* 



x N 



let L denote the autocorrelation function of K, videlicet 

oo 

nL{x) = J (I + f)- 2 cos{x£) d£ 



for — oo < x < oo. We have 



G 



L 



Xj Xfc 



_£(14) / x 3 ~ x k 



j,k=l,...,N 



j,k=l,...,N 

L(x) = ? - K(x) - I K'(x), L^\x) = I K^(x) + | K^(x) 
for — oo < x < oo, and 

(C.8) S.R.Q. = least eigenvalue of H with respect to G 

— in other words, 

S.R.Q. = least eigenvalue of G'^H G' 1 ^ . 
Figure 12 shows plots of S.R.Q. versus (6 — a)/((N — The following table 



N 


minimum value 


(6- 


-a)/((N-l)[x) 




[i/(b - a) 


5 





115320757000 





900127730000 





277738360532 


10 





042825969700 





578865888000 





191946206219 


20 





015529762000 





361101444000 





145752889726 


30 





008487160550 





271108504000 





127191726235 


40 





005502126710 





220171280000 





116459447577 


50 





003920874640 





187117585000 





109065982576 


60 





002967795760 





163520070000 





103651818045 


70 





002342495540 





145755736000 





099431789245 


80 





001906845850 





131825897000 





096022315313 


100 





001349762400 





111466361000 





090619358257 


120 





001016274330 





097183852500 





086468699567 


140 





000798662701 





086282601400 





083380015062 


150 





000716797334 





081805225800 





082041328416 
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lists sample solutions to equation (|C.7j) . 
The following formula 

(C.9) — — w 0.0415 + 0.5416 x AT 1 / 2 _ 0.6426 x A^ 1 + 1.3706 x A^ 3 / 2 
o — a 

gives an effectual estimate of such solutions. 
C.4. The following theorem holds. 

Theorem C.l. Suppose u solves problem (|C.1|) & (|C.2[) ; suppose f,e and E obey 

\f(x k )- 9k \<e (k = l,...,N), 



M J fix) >dx + fi J [f""(x)Ydx < 2(6 - a)~ 6 E z ; 

— oo — OO 

let 

* = max{(l-l/A0- 1 / a J s (JV-l)- 2 }. 
If S approaches zero, and 

A / fi \ 3 / £ N ~ 2 



N \b~aj \E, 
is bounded and bounded away from zero, then 

E- 1 max{|/0) - u(a;)| : a < x < b} = C>(<5 3/4 ), 

ET 1 ^ - a) max{|/'(x) - u'(x)\ : a < x < b} = 0(S 1/4 ). 

Proof. Let 

(CIO) v = f-u. 

Functional J, whose domain is W A,2 (—oo, oo) and whose value at any trial func- 
tion ip obeys 

N °° 

%) = ^[^)-rf+A / [^ 7 (0 2 + M-V 2 ] dx, 

k— 1 

attains its minimum value at u. Consequently, 

J(f) = J{u) + (a remainder), 

N °° 

remainder = Y^H^k)] 2 + A / [pJ {v""f + pL^v 2 } dx. 



Since 



k=l 

J(f)<Ns 2 + 2X \ 2 , 
J(u) > 0, 

OO OO 

J [m 7 (u"") 2 + M" 1 " 2 ] dx>2(i 3 J {v") 2 dx, 



we infer 



n y f 2 

(C.ll) J2[v(xk)} 2 + 2\^ (v") 2 dx<E 2 lN^ +2A 



\ 3' 



b — a . 
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Combining inequality (jC.llj) and Lemma IC.2I below results in 
b 

(b-a)- 1 I v 2 dx < 
E 2 



N 1 N ( (i Y 3 \ j / e\ 2 2A / n x ' 



N — X tt 4 {N - l) 4 2A \b-a 
inequality (|C.11|) also yields 



E) N \b — a 



[b-af J{v"fdx<E 2 jl + g 

n > 



N ( n \ " t e \ 2 



6-a/ \E 



The last two inequalities and a hypothesis imply 



(C.12) 
(C.13) 



(b -a)- 1 J v 2 dx < (Const.) E 2 S 2 

a 

{b-af}{v") 2 dx < (Const.) E 2 



The conclusions follow from (|C.10p . (IC.12|1 and (|C.13|) . by virtue of Lemma IC31 
below. □ 




J o ° o 



0.1 0.2 



0.4 0.5 0.6 0.7 0.8 0.9 



Figure 13. Original function (dotted), and samples polluted by 
noise (circled). 



Lemma C.2. Let —oo < a < b < oo, and N — 2,3, ... ; let 
b — a 



Ax 



N - V 



x k =a+(k-l)Ax (k = l,...,N). 



:!') 



The following inequality 

1/2 



Jv 2 dx\ <{Ax)- l '*^pv(x k ) a *cj U(v") 2 

holds for every v from W 2,2 (a, b). 

Lemma C.3. Suppose v is in W 2,2 (a, b), and 

(b - a)- 1/2 \ [ v 2 dx } =\\v\\, (b-af f(v") 2 dx = 1. 



1/2 



The following inequalities hold 

max M < 2 x /4 3- 3 / 8 ||w|| 3 / 4 + 0(||»||), 



{b -a) max\v'\ < 2 1 / 4 3- 3 /*\\v\\V* + 0(\\v\\). 




0.1 0.2 0.3 0.4 0.5 



0.7 0.8 0.9 



FIGURE 14. Original and recovered functions. 



The proofs of Lemma IC.2I and IC.3I are beyond the scope of the present paper, 
and are omitted. 

C.5. Here is an example, demonstrating how the present method works. Let 

a = 0, 6=1, 
f{x) = exp [-72 (x - I) 2 )} [cos(25a;) - 4 sin(25x)], 
N = 150, 

9k = f( x k) ± (5% random noise) (k — 1, .... N). 
Let parameter A obeys 

A= 1(T 6 , 
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80 




qI 1 1 1 1 1 1 1 1 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



Figure 15. Original and recovered derivatives. 

let parameter [i be given by formula (|C.9|) . and let u solve problem (|C.1|) & (|C.2[) . 
Figure 13 plots g±,g2, ■ ■ ■ ,9n versus Xi,x%, ■ ■ ■ ,xn\ figure 14 plots / and u, and 
figure 15 plots /' and u' . 
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